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I. INTRODUCTION 



The last two decades have been witness to intense activities in experimental and theo- 
retical study of correlation properties of interacting quantum systems. The problem is still 
challenging and open in the quantum many-body systems which are non-integrable. The task 
becomes simplified for integrable systems like one dimensional uniform Bose gas described 
by the Lieb-Liniger (LL) model [1,2], which assumes that particles interact via a 5-function 
repulsive potential. Since the experimental observation of Bose Einstein condensate (BEC) 
in ultracold trapped alkali atomic vapors, a lot of theoretical and experimental work has 
been done to study its correlation properties [3-12]. In the recent experimental situation, 
it is easy to achieve a quasi-one- dimensional strongly interacting degenerate Bose gas in a 
highly anisotropic trap [13-15] which is correctly described by the Lieb-Liniger model [16]. 
The complementary case is the weakly interacting (via finite range forces) trapped Bose 
gas. Although it is commonly believed that the Gross Pitaveskii (GP) equation is adequate 
for weakly interacting Bose gases, but a more rigorous and accurate many-body treatment, 
incorporating realistic interatomic interactions and interatomic correlations, is crucial for 
studying correlation properties in realistic condensates. The first motivation of the present 
work is to investigate the importance of finite size and trapping effects on the ground state 
correlation properties of ultracold atomic BECs. The second motivation of our study is to 
assess the validity of shape-independent approximation in the correlation properties. Dilute 
BECs are known to possess shape-independent (SI) property which is frequently described 
by the 5-function potential, whose strength is proportional to the s-wave scattering length 
(a s ). In the mean-field description, the effective interaction is determined by the factor 
Na s (where a s is expressed in units of oscillator length of the trap), whereas in the full 
many-body description we solve the many-body Schrodinger equation for different two-body 
potentials that generate identical a s . Thus, our present calculation serves as a stringent 
test of the SI approximation and verifies whether the long range correlation at all affects 
the correlation properties or not. Inclusion of all possible two-body correlations in our ab 
initio many-body method correctly calculates correlation effects at zero temperature. Use 
of the realistic van der Waals interatomic potential with a long attractive tail provides the 
realistic aspects of the correlation function and its momentum distribution. Specially the 
effect of short range repulsion in the interatomic potential is expected to be reflected in the 
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pair correlation function. 



The paper is organized as follows. In the next section, we briefly review our theoretical 
approach. In Section III, we present the results of our calculation of one-body density, pair 
distribution function and their momentum distributions. Finally, in Section IV we draw our 
conclusions. 

II. POTENTIAL HARMONIC EXPANSION METHOD 

We adopt the potential harmonics expansion method together with a short range cor- 
relation function (CPHEM), which has already been established as a successful and useful 
technique for investigating BEC with realistic two-body interactions. In the following, we 
briefly describe the technique. Interested readers can find the details in Refs.— ~— . 

We start with the many-body Schrodinger equation of A spinless bosons trapped in an 
isotropic harmonic oscillator potential of frequency u at zero temperature. The center of 
mass motion can be separated by introducing the center of mass vector (R) and N = A — 1 
Jacobi vectors ...,(n} defined as^ 



where Xi is the position vector of the z-th particle. The relative motion of the system is 
described by 



where V trap and V int are respectively the trapping and pair-wise interaction potentials, ex- 
pressed in terms of the Jacobi vectors. The energy of the relative motion is Er. Next 
hyperspherical variables are introduced by defining a 'hyperradius' 
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and a set of (3N — 1) 'hyperangles', consisting of 2N polar angles of N Jacobi vectors and 
(N — 1) angles defining their relative lengths^. In the hyperspherical harmonics expansion 
method (HHEM) the iV-body Schrodinger Equation ([2]) is solved by expanding ip in the 
complete set of hyperspherical harmonics (HH), which are the eigenfunctions of the N- 
dimensional Laplace operator (analogous to the spherical harmonics in three dimension)^. 
Substitution of this in Eq. ([2]) and projection on a particular HH results in a set of coupled 
differential equations. However imposition of symmetry of the wave function and calculation 
of the matrix elements are such formidable tasks that a practical solution for A > 3 is nearly 
impossible. Moreover due to the fact that the degeneracy of the HH basis increases very 
rapidly with the increase in the grand orbital quantum number K-, a convergent calculation 
using HHEM with the full HH basis is extremely computer intensive. This is the price one 
pays for keeping all many-body correlations. On the other hand, a typical experimental 
BEC is designed to be extremely dilute to eliminate the possibility of molecule formation 
through three-body collisions. Hence we assume that the probability of three and more 
particles to come within the range of interatomic interaction is negligible and the effect of 
the two-body correlations will be adequate for the full many-body wave function. Moreover, 
only two-body interactions are relevant. Therefore, ip can be expressed in terms of Faddeev 
components ipij of the (ij) interacting pair- 



Since only two-body correlations are important, ipy is a function of the interacting-pair sep- 
aration Tij = Xi — Xj and r only. Hence it can be expanded in the subset of HH needed 
for the expansion of the two-body interaction Vffij). This subset is called the potential 
harmonics (PH) basis^. This results in a dramatic simplification in the analytic and com- 
putation works and it has been used upto 14000 bosons in the trap>^. In this procedure, a 
realistic interatomic interaction can be used. Such an interaction has a strongly repulsive 
core. Hence ipij must be vanishingly small for small values of r^. But the first (K = 0) 
term of the PH basis is a constant^. Hence to reproduce the short-range behaviour of ifi^ 
correctly, a large number of PH basis functions is needed, which slows down the rate of 
convergence considerably. Hence, we introduce a 'short-range correlation function' in the 
PH expansion basis^. In the zero temperature BEC, the kinetic energy of the interacting 
pair is practically zero. Therefore, the two-body collision process is determined entirely by 
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the s-wave scattering length (a s ). The small behavior of ipij will be that of the zero- 
energy wave function rj(rij) of the pair interacting via V(rij). Hence we introduce rj(rij) as 
a short-range correlation function in the PH expansion^. We have checked explicitly that 
this improves the rate of convergence dramatically. In order that the small behavior of 
4 > ij{ r ij-i r ) corresponds to the correct two-body interaction appropriate for the experimental 
a s , we adjust the very short-range behavior of the interatomic interaction (in the case of 
the van der Waals potential used by us, it will be the hard-core radius r c ) so that rj(rij) 

asymptotically becomes proportional to (1 Substitution of this expansion of ip in 

Eq. ([2]) and projection on a particular PH gives rise to a set of coupled differential equations 
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in v ±s. 



where £ = Z + (3A-6)/2, U Kl (r) = f K iu l K (r), a = and = l + Here, / is the orbital 
angular momentum of the system (assumed to be contributed by the interacting pair only) 
and u l K (r) is the coefficient of expansion of ipij(rij,r) in the correlated PH basis, while fxi 
represents the overlap of the PH corresponding to the (ij) -partition with the sum of PHs of 
all partitions^. Expressions for the potential matrix element VKK>{ r ) and fxi can be found 



Eq. (5) is solved by hyperspherical adiabatic approximation (HAA)^, for which we assume 
that the hyperangular motion is much faster than the hyperradial motion. Consequently, 
the former can be separated adiabatically and solved to obtain an effective potential as a 
parametric function of r, in which the hyperradial motion takes place. The hyperangular 
motion is effectively solved by diagonalizing the potential matrix Vkk 1 if) together with the 
hyper-centrifugal potential of Eq. (5) to get the lowest eigenvalue wo(r) [corresponding eigen 
column vector being Xxo( r )L which is the effective potential for the hyperradial motion. 
Finally, the adiabatically separated hyperradial equation 
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is solved subject to appropriate boundary conditions to get Er and the hyperradial wave 
function Co( r )- The many-body wave function can be constructed in terms of Co( r ) and 



III. RESULTS 



Now in order to study the correlation properties of the interacting Bose gas we choose a 
realistic interatomic potential having a strong repulsive core at a small separation and an 
attractive tail at large atomic separations. This is approximately represented by the van 
der Waals potential with a hard core of radius r c and a — attractive tail, viz. V(rij) = oo, 

for Tij < r c and — =- for > r c . The effective interaction is characterized by the s-wave 



scattering length (a s ), which depends strongly on r^. Usually the potentials are chosen 
to be purely attractive or purely repulsive according to whether a s is negative or positive 
respectively. In our many-body calculation, we solve the zero-energy two-body Schrodinger 
equation with Vfcj) to obtain aj^^i. The value of r c is adjusted so that a s has the values 
corresponding to the JILA experiments^ 1 ^. A typical value of r c is of the order of 10 -3 
o.u. In atomic units this is a few tens of Bohr, which is larger than the atomic radius. Note 
that r c is expected to be larger than the atomic radius, as the van der Waals potential with 
the hard core of radius r c is the effective potential which produces the correct experimental 
zero- energy scattering cross- section (given by a s ). 

A. One-body density 

We define the one-body density, as the probability density of finding a particle at a 
distance from the center of mass of the condensate 



where ip is the full many-body wave function and the integral over the hypervolume r' 
excludes the variable r&. After a lengthy but straightforward calculation we arrive at a 
closed form given by 
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where D = 3 A — 3 and /i^ is the norm of the Jacobi polynomial P^ p (z). 



a/3. 



The one-body density contains information regarding one particle aspect of the bosonic 
system. Although it is not directly measurable but in the interferometry experiment one 
can indirectly explore itr^. In Fig.l we present calculated one-body density as a function of 
the distance from the trap centre for a repulsive interaction corresponding to a s = .00433 
o.u. for 100 87 Rb atoms in the JILA trap^. In our calculation, length and energy are 
measured in oscillator units (o.u.) of length (aho = y^jj an d energy {ho) respectively. 
For comparison, we also include the mean-field results. The effect of interaction is revealed 
by the deviation from the Gaussian profile. To explore the effect of interaction, we also 
calculate one-body density for a smaller a s value viz., a s = 2.09 x 10 -4 o.u. for the same 
number of particles. Since the trapped condensate is always stable even for a large A, we 
see appreciable changes in Ri(ff.) as a s decreases. For small a s , the density distribution is 
sharper as the correlations induced by the interactions are weak, while for the larger a s , the 
peak is flatter with a larger width. In Fig. 2 we present calculated one-body density for 
an attractive BEC with a s = —1.836 x 10~ 4 o.u for different number of 85 Rb atoms in the 
JILA trap2& The effective interaction parameter is A = A\a s \/a>h . For A— 100, A is small 
and the system exhibits weak one-body density which extends to the trap size. Increasing 
A gradually, the net effective attraction becomes strong and the density becomes sharply 
peaked at a smaller distance. 

Taking the Fourier transform of Ri(fk), we obtain the one-body momentum distribution. 
One-body density is an abstract concept. However, its fourier transform gives the experi- 
mentally measurable quantity, the momentum distribution. In Fig. 3, we plot our results 
for repulsive interactions. It evolves from a Gaussian in the noninteracting limit to a curve 
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FIG. 1. (Color online) One-body density distribution as a function of (in o.u.) for a repulsive 
BEC with A=10000 bosons. The choice of a s = 0.00433 o.u. corresponds to 87 Rb experiment 
in the JILA trap. PHEM corresponds to our present many-body results and GP corresponds to 



mean-field results. 
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FIG. 2. (Color online) One-body density distribution as a function of (in o.u.) for an attractive 
interaction (a s = —1.836 x 10 -4 o.u. for 85 Rb atoms in the JILA trap), for various indicated values 
of particle numbers. 



having a sharper peak, as the net interaction increases. The peak at k = becomes more 
pronounced with increase in effective repulsion, whereas for weak interaction it develops a 
long-range tail in the momentum space. The momentum is being redistributed to higher k 
values. The width of the low-momentum peak for a s = 0.00433 o.u. and A = 10000 is about 
0.7 /mi -1 . 



We have remarked earlier that with recent progress in creating atomic clouds with large 
dipole moment, interest has been shifted to longer range interaction instead of taking only 
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FIG. 3. (Color online) Calculated one-body momentum distribution for the repulsive BEC for 
different values of A and a s . The choice of a s = 0.00433 o.u. corresponds to 87 Rb experiment in 
JILA trap. All quantities are in appropriate oscillator units. 

contact interaction. Thus, even in the low-density limit the use of a realistic interatomic in- 
teraction potential in the many-body calculation has been emphasized by several authors [27- 
29] . The van der Waals potential is an ideal choice as it properly takes care of the effect of 
realistic dipole-dipole interaction. The strength of the van der Waals interaction C§ changes 
widely from a small value for H atoms to a high value for Cs atoms. Hence, to see how the 
one-body density is affected by the strength of the long-range tail of the two-body potential, 
we select three different values of C§ in addition to the actual experimental value - one below 
and two above. These are (in o.u.): 5 x 10 -11 , 6.489755 x 10" 11 , 8.5 x 1CT 11 , 8 x 1CT 10 , the 
second one being the experimental value. For each value of Cq, we calculate corresponding 
r c as before, such that a s has the experimental value 0.00433 o.u. For each set of (C$, r c ), 
we solve the many-body Schrodinger equation as before. Calculated one-body density for 
different sets are shown in Fig. 4 for 10000 bosons. We observe that the one-body den- 
sity of the inhomogeneous gas for different two-body potentials are almost indistinguishble; 
i.e. independent of the shape of V(ry). Our numerical calculations confirm that one-body 
density is absolutely determined (within numerical errors) by the parameter a s only. 

This is in contrast with ref.— , where ground state energy of the condensate was found to 
depend on the shape of the potential. Thus the one-body density does not depend strongly 
on the shape of the interaction potential, while the ground state energy does. 

An explanation of the above observation is as follows. The total condensate energy 
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FIG. 4. (Color online) One-body density distribution for the repulsive BEC with several Cq pa- 
rameters, all of which correspond to the identical a s = 0.00433 o.u., with 10000 87 Rb atoms in the 
JILA trap. All quantities are in oscillator units. 

depends on the minimum of the well, in addition to the actual shape of the lowest eigen 
potential. In Table 1, we present the dependence of the position (r m ) and the value (co 0m ) of 
the minimum of lowest eigen potential for different values of C§ parameter, corresponding 
to the same a s . 

Table 1. Dependence of r m and u 0m on C 6 (corresponding to a s = 0.00433) for A = 10000 
atoms (all quantities are expressed in o.u.). 
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We observe that, as C 6 increases, u 0m decreases appreciably, while position of the minimum 
of effective many-body potential changes by a small amount. In the next sub-section, we 
will see that the shape of the ooq(v) curve remains practically unchanged. As the total 
condensate energy depends on the depth of the potential, as well as on its shape (stiffness), 
we observe that the ground state energy decreases gradually with increase in C^L. Thus 
the shape independence hypothesis is violated for the total energy. However, the one-body 
density distribution is given by the many-body wave function, which is independent of the 
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minimum of the potential, but depends only on its shape. Hence it remains unchanged with 
change in Cq parameter and shape independence of the one-body density profile is satisfied. 

B. Pair distribution function 

Another key quantity is the pair distribution function R^rij), which determines the 
probability of finding the (zj)-pair of particles at a relative separation r^. The study of 
pair correlation is important, since the interatomic interactions play a crucial role as there 
are two competing interaction length scales. When the atoms try to form clusters, a strong 
very short range repulsion in atomic interaction comes into play. As atoms repel each other 
strongly at very small separations it is impossible to get some non- vanishing value of ^(fy) 
at Tij = 0. We calculate it as 

Mru) = [ \^\ 2 dr", (9) 

where i/j is the many-body wavefunction and the integral over the hypervolume r" excludes 
integration over r^. Again after a lengthy calculation we can put it in a closed form given 
by 




XKo(r) X K'o(r)Pf(z)P%(z)dz. (10) 

We plot R,2{rij) in Fig. 5 for the attractive interaction with different particle numbers. 
As mentioned earlier, pair correlation vanishes cLS T*2 j r 0, due to the strong interatomic 
repulsion and it cannot extend beyond the size of the condensate. Hence R2 is peaked 
at some intermediate value of r^. For weak interactions (small A), the correlation length 
is large. However as the effective attractive interaction increases, two interacting particles 
come closer and the pair-correlation length decreases; the pair distribution function becomes 
sharply peaked. It indicates stronger pair correlation in the system, in agreement with 
expectations. However our results are at variance with those obtained in LL model^, which 
describes uniform systems with no confinement. Consequently pair correlation function 
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approaches its maximum asymptotically, whereas in our confined three dimensional case, it 
vanishes asymptotically. When A is increased to 2470, which is very close to the critical 
number ( A cr ~ 2475), the effect is quite prominant. As the width of the curves decrease 
with increase in effective attractive interaction, the interacting pair becomes more localized. 
This manifests the possiblity of clustering due to large two-body interaction and three- 
body recombination, if the attraction becomes strong enough. This will lead to an eventual 
collapse of the attractive condensate, when the pair-correlation length will drastically reduce 
to a very small value and the condensate will be destroyed. 
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FIG. 5. (Color online) Plot of pair distribution function (i?2( r ij)) against rjj (in o.u.) for an 
attractive Bose gas with different particle numbers in the trap. 

Two-particle correlation is directly related with pair inelastic processes and can be di- 
rectly measured in photoassociation in interatomic collisions. As the attractive interaction 
increases, the gas becomes highly correlated and different inelastic processes may take place. 
The opposite case is the repulsive Bose gas where the pair inelastic process will be suppressed 
with increase in interaction. For completeness, in Fig. 6 we plot the relative momentum 
distribution obtained as the Fourier transform of i?2(^) for various particle numbers for 
87 Rb condensate in the JILA trap, with a s = 0.00433 o.u. (100 a.u.). The width of the 
curve decreases and exhibits a sharp peak at k = 0, whose height increases with increasing 
number of particles. 

Lastly to visualize the effect of different strength parameter C§ on the pair-correlation, 
we plot -R2( r ij) against for the previous sets of potential parameters in Fig. 7. We find 
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FIG. 6. (Color online) Plot of relative momentum distribution against k (in (o.u.) -1 )for vari- 
ous indicated values of particle numbers. The choice of a s = 0.00433 o.u. corresponds to 87 Rb 
experiment in JILA trap. 



that different interatomic potential shapes corresponding to the same a s produce identical 
correlation properties, which indicates that the calculated property is independent of the 
shape of the two-body potential. The fact that the ground state energy is dependent on the 
shape of the two-body potential whereas the one-body density and the pair distribution 
function are not, can be attributed to the fact that the many-body wave function remains in- 
variant irrespective of the variation of the shape of the two-body potential due to variation of 
(Cq, r c ) for the same value of a s . Note that the condensate wave function in the hyperradial 
space Co( r ) is obtained by solving the hyperradial equation (6) in the effective many-body 
potential 0Jo(r). The shape and stiffness of this effective potential remain unchanged for the 
same value of a s , irrespective of the variation of the shape of the two-body potential. How- 
ever, as discussed earlier, the position and value of the minimum of the effective many-body 
potential change with the variation of the parameter Cq. To demonstrate this, we plot the 
effective many-body potential 0Jo(r) as a function of r for various values of Cq in Fig. 8. The 
curves have been shifted both vertically and horizontally in order that the minima coincide 
(both in position and value). One observes that all the curves overlap completely (within 
numerical errors) showing that the shape and stiffness of the effective potential remain in- 
variant. Hence the kinetic energy and the wave function remain unchanged, while potential 
energy and the total energy change with change of Cq. By direct numerical calculation, we 
have checked that the change of kinetic energy is very small compared with that of the total 
energy, as Cq varies. Therefore the ground state energy is shape dependent, whereas the 
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one-body density and pair distribution function are not. 




FIG. 7. (Color online) Plot of two-body correlation function against rij in o.u. for repulsive BEC 
with several indicated Cq parameters, all of which correspond to identical a s = 0.00433 o.u. with 
10000 87 Rb atoms in JILA trap. (The portion near the origin is magnified in the inset to show 
that the pair distribution function vanishes inside the hard core. The same is true for Fig. 5.) 



IV. CONCLUSION 

In conclusion, the present work focuses on the one-body density and pair-correlation 
aspects of the zero temperature weakly interacting trapped Bose gas. The use of correlated 
many-body approach takes care of the effect of finite-size, where quantum fluctuation is 
important and gives a realistic picture of correlation properties. Due to the use of a realistic 
interatomic interaction and consideration of finite number of atoms in the trap, our results 
deviate from the earlier results^, but exhibit realistic aspects which are relevant to exper- 
iments. Our calculation is performed for a two-body potential (van der Waals potential), 
whose parameters can be adjusted to give both positive and negative scattering lengths, for 
Rb atoms in the JILA trap. Thus, our results are realistic and can be experimentally verified 
in future. Our calculations also verify the validity of the shape-independent approximation 
in dilute BECs for one-body density and pair-distribution functions. This is in contrast 
with the earlier observed shape dependence for ground state energies^. 
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FIG. 8. (Color online) Plot of the effective potential (in o.u.) against r (in o.u.) for 10000 87 Rb 
atoms with a s = 0.00433 o.u.. The curve corresponding to Cq = 5.0 x 10~ n is in its actual position. 
The other curves have been shifted both vertically and horizontally to coincide the minima (both 
in position and value). All the curves overlap completely (within numerical errors). This shows 
that the shape of the effective many-body potential u>o(r) is independent of the choice of (Ce,r c ) 
for the same value of a s . 
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